cap program drop distance

program define distance
syntax varlist (min=4 max=4 numeric) [if] [in] [, Metric(string) Aname(string) Bname(string) DIRection TN]
// IMPORTANT: 1: latitude A, 2: longitude A, 3: latitude B, 4: longitude B


// parameters
local conv 	= _pi / 180

// metric
if inlist("`metric'","mi") {
	local radius = 6378.388 * 0.621371192
	local metstr "miles"
}
else {
	local radius = 6378.388
	local metstr "kilometers"
}


//names
if "`aname'"=="" 	local Al "A"
else 				local Al "`aname'"
if "`bname'"=="" 	local Bl "B"
else 				local Bl "`bname'"

local A = upper(substr("`Al'",1,2))
local B = upper(substr("`Bl'",1,2))


// Quelle: http://de.wikipedia.org/wiki/Orthodrome#Winkelberechnung

// latitude/longitude of A and B in radian
gen lat_`A' 	= `1' * `conv'
gen long_`A' 	= `2' * `conv'
gen lat_`B' 	= `3' * `conv'
gen long_`B' 	= `4' * `conv'

// distance in radian
gen double dist = acos(sin(lat_`A') * sin(lat_`B') + cos(lat_`A') * cos(lat_`B') * cos(abs(long_`A' - long_`B')))

// distance in km or mi
gen double dist_`A'`B' = (`radius') * dist
lab var dist_`A'`B' "Distance: `Al'-`Bl' (in `metstr')"

// angle of direction
if "`direction'"!="" {
	gen double dir_`A' 	= (1/`conv') * acos((sin(lat_`B') - sin(lat_`A')*cos(dist)) / (cos(lat_`A')*sin(dist)))
	gen double dir_`B' 	= (1/`conv') * acos((sin(lat_`A') - sin(lat_`B')*cos(dist)) / (cos(lat_`B')*sin(dist)))
		
	lab var dir_`A' "Angle in `A'"
	lab var dir_`B' "Angle in `B'"	
	
	// true north angle
	if "`tn'"!="" {	
	
		gen double 	tn_`A'`B' = dir_`A' 		if long_`A' < long_`B'
		replace    	tn_`A'`B' = 360 - dir_`A' 	if long_`A' > long_`B'	
		
		gen double 	tn_`B'`A' = dir_`B' 		if long_`B' < long_`A'
		replace 	tn_`B'`A' = 360 - dir_`B'	if long_`B' > long_`A'		
		
		lab var tn_`A'`B' "True north: `Al' --> `Bl'"
		lab var tn_`B'`A' "True north: `Bl' --> `Al'"
		
		*drop dir_`A' dir_`B'
	}
}



drop dist lat_`A' long_`A' lat_`B' long_`B'

end
